arXiv:1505.01429v3 [cs.CV] 5 Jan 2016 


1 


Geometry-Aware Neighborhood Search for Learning 
Local Models for Image Reconstruction 
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Abstract —Local learning of sparse image models has proven 
to be very effective to solve inverse problems in many computer 
vision applications. To learn such models, the data samples 
are often clustered using the K-means algorithm with the Eu¬ 
clidean distance as a dissimilarity metric. However, the Euclidean 
distance may not always be a good dissimilarity measure for 
comparing data samples lying on a manifold. In this paper, 
we propose two algorithms for determining a local subset of 
training samples from which a good local model can be computed 
for reconstructing a given input test sample, where we take 
into account the underlying geometry of the data. The first 
algorithm, called Adaptive Geometry-driven Nearest Neighbor 
search (AGNN), is an adaptive scheme which can be seen as an 
out-of-sample extension of the replicator graph clustering method 
for local model learning. The second method, called Geometry- 
driven Overlapping Clusters (GOC), is a less complex nonadap- 
tive alternative for training subset selection. The proposed AGNN 
and GOC methods are evaluated in image super-resolution, 
deblurring and denoising applications and shown to outperform 
spectral clustering, soft clustering, and geodesic distance based 
subset selection in most settings. 

Index Terms —Clustering, patch manifolds, nearest neighbor 
search, image superresolution, image deblurring, image restora¬ 
tion. 


1. Introduction 

M any image restoration problems such as superresolu¬ 
tion, deblurring, and denoising can be formulated as 
a linear inverse problem, by modeling the image deformation 
via a linear system. Such problems are generally ill-posed and 
the solutions often rely on some a priori information about 
the image to be reconstructed. Research in the recent years 
has proven that adopting an appropriate sparse image model 
can yield quite satisfactory reconstruction qualities. Sparse 
representations are now used to solve inverse problems in 
many computer vision applications, such as superresolution 
a PL 0; denoising Q compressive sensing 

0’ 13’ deblurring |[Tj7 ||2]f While several works 

assume that the image to be reconstructed has a sparse 
representation in a large overcomplete dictionary 0, 0, it has 
also been observed that representing the data with small, local 
models (such as subspaces) might have benefits over a single 
and global model since local models may be more adaptive 
and capture better the local variations in data characteristics 
0, 0, fTQ| . The image restoration methods in 0 and 0 
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propose a patch-based processing of images, where the training 
patches are first clustered and then a principal component 
analysis (PCA) basis is learned in each cluster. The idea of 
learning adaptive models from groups of similar patches for 
image restoration has been exploited in several recent works 
(TT), (T2), (T3) 

When learning local models, the assessment of the similarity 
between image patches is of essential importance. Different 
similarity measures lead to different partitionings of data, 
which may eventually change the learned models significantly. 
Many algorithms constructing local models assess similarity 
based on the Euclidean distance between samples. For example 
in 0 and 0 image patches are clustered using the K-means 
algorithm, where patches having a small Euclidean distance 
are grouped together to leam a PCA basis. Test patches are 
then reconstructed under the assumption that they are sparsely 
representable in this basis. 

However, patches sampled from natural images are highly 
structured and constitute a low-dimensional subset of the high¬ 
dimensional ambient space. In fact, natural image patches are 
commonly assumed to lie close to a low-dimensional manifold 
|T^, |T5| . Similarly, in the deconvolution method proposed 
in |10|, image patches are assumed to lie on a large patch 
manifold, which is decomposed into a collection of locally 
linear models learned by clustering and computing local PCA 
bases. The geometric structure of a patch manifold depends 
very much on the characteristics of the patches constituting 
it; the manifold is quite nonlinear especially in regions where 
patches have a rich texture. When evaluating the similarity 
between patches on a patch manifold, care should be taken 
especially in high-curvature regions, where Euclidean distance 
loses its reliability as a dissimilarity measure. In other words, 
in the K-means based setting of 0 and 0, one may obtain a 
good performance only if the local PCA basis agrees with the 
local geometry of the patch manifold, i.e., the most significant 
principal directions should correspond to the tangent directions 
on the patch manifold so that data can be well approximated 
with a sparse linear combination of only a few basis vec¬ 
tors. While this easily holds in low-curvature regions of the 
manifold where the manifold is fiat, in high-curvature regions, 
the subspace spanned by the most significant principal vectors 
computed from the nearest Euclidean-distance neighbors of 
a reference point may diverge significantly from the tangent 
space of the manifold if the neighborhood size is not selected 
properly p^ , (Tv). This is illustrated in Figure 0 where the 
first few significant principal directions fail to approximate 
the tangent space because the manifold bends over itself as 
in Figure |l(b)| or because the curvature principal components 
dominate the tangential principal components as in Figure [T^ 
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In this work, we focus on image restoration algorithms 
solving inverse problems based on sparse representations of 
images in locally learned subspaces, and we present geometry- 
driven strategies to select subsets of data samples for learning 
local models. Given a test sample, we address the problem 
of determining a local subset of the training samples, i.e., a 
neighborhood of the test sample, from which a good local 
model can be computed for reconstructing the test sample, 
where we take into account the underlying geometry of the 
data. Hence, the idea underlying this work is to compute local 
models that agree with the low-dimensional intrinsic geometry 
of data. Low dimensionality allows sparse representations of 
data, and the knowledge of sparsity can be efficiently used for 
solving inverse problems in image restoration. 

Training subsets for learning local models can be deter¬ 
mined in two ways; adaptively or nonadaptively. In adaptive 
neighborhood selection, a new subset is formed on the fiy 
for each test sample, whereas in nonadaptive neighborhood 
selection one subset is chosen for each test sample among 
a collection of training subsets determined beforehand in 
a learning phase. Adaptive selection has the advantage of 
fiexibility, as the subset formed for a particular test sample 
fits its characteristics better than a predetermined subset, but 
the drawback is the higher complexity. In this work, we study 
both the adaptive and the nonadaptive settings and propose two 
different algorithms for geometry-aware local neighborhood 
selection. 

We first present an adaptive scheme, called Adaptive 
Geometry-driven Nearest Neighbor search (AGNN). Our 
method is inspired by the Replicator Graph Clustering (RGC) 
p8| algorithm and can be regarded as an out-of-sample exten¬ 
sion of RGC for local model learning. Given a test sample, the 
AGNN method computes a diffused affinity measure between 
each test sample and the training samples in a manner that 
is coherent with the overall topology of the data graph. The 
nearest neighbor set is then formed by selecting the training 
samples that have the highest diffused affinities with the test 
sample. 

The evaluation of the adaptive AGNN method in super¬ 
resolution experiments shows a quite satisfactory image re¬ 
construction quality. We then propose a nonadaptive scheme 
called Geometry-driven Overlapping Clusters (GOC), which 
seeks a less complex alternative for training subset selection. 
The method computes a collection of training subsets in a 
prior learning phase in the form of overlapping clusters. The 
overlapping clusters are formed by first initializing the cluster 
centers and then expanding each cluster around its central 
sample by following the AT-nearest neighborhood connections 
on the data graph. What really determines the performance 
of the GOC method is the structure of the clusters, driven by 
the number of neighbors K and the amount of expansion. We 
propose a geometry-based strategy to set these parameters, by 
studying the rate of decay of PCA coefficients of data samples 
in the cluster, thereby characterizing how close the cluster lies 
to a low-dimensional subspace. 

Note that, while the proposed AGNN and GOC algorithms 
employ similar ideas to those in manifold clustering methods. 


our study differs from manifold clustering as we do not aim 
to obtain a partitioning of data. Instead, given a test sample to 
be reconstructed, we focus on the selection of a local subset 
of training data to learn a good local model. We evaluate 
the performance of our methods in image superresolution, 
deblurring and denoising applications. The results show that 
the proposed similarity assessment strategies can provide 
performance gains compared to the Euclidean distance, espe¬ 
cially for superresolving images with rich texture where patch 
manifolds are highly nonlinear. When applying the proposed 
method in the superresolution problem, we select the NCSR 
algorithm Q as a reference method, which currently leads 
the state of the art in superresolution. We first show that 
the proposed AGNN and GOC methods outperform reference 
subset selection strategies such as spectral clustering, soft clus¬ 
tering, and geodesic distance based neighborhood selection. 
Finally, we perform comparative experiments with the NCSR 
Q, ASDS Q, and SPSR |T9| superresolution algorithms, 
which suggest that the proposed methods can be successfully 
applied in superresolution for taking the state of the art one 
step further. The experiments on image deblurring also confirm 
these findings, suggesting that the proposed methods perform 
better than K-means in most images. Meanwhile, we have 
achieved a marginal performance gain in image denoising 
applications only at small noise levels. 

The rest of the paper is organized as follows. In Section |n| 
we give an overview of manifold-based clustering methods. In 
Section III we formulate the neighborhood selection problem 
studied in this paper. In Section [Tvl w e discuss the proposed 
AGNN method. Then in Section M we describe the GOC 
algorithm. In Section VI we present experimental results, and 
in Section (Vnl we conclude. 


II. Clustering on manifolds: related work 

As our study has close links with the clustering of low¬ 
dimensional data, we now give a brief overview of some 
clustering methods for data on manifolds. The RGC method 
|T^ , from which the proposed AGNN method has been 
inspired, first constructs a data graph. An initial affinity matrix 
is then computed based on the pairwise similarities between 
data samples. The affinity matrix is iteratively updated such 
that the affinities between all sample pairs converge to the 
collective affinities that consider all paths on the data graph. 
Spectral clustering is another well-known algorithm for graph- 
based clustering p0| , Samples are clustered with respect 
to a low-dimensional embedding given by the functions of 
slowest variation on the data graph, which encourages assign¬ 
ing neighboring samples with strong edge weights to the same 
cluster. The Laplacian eigenmaps method p2| builds on the 
same principle; however, it targets dimensionality reduction. 

Geodesic clustering provides an extension of the AT-means 
algorithm to cluster data lying on a manifold, where the 
Euclidean distance is replaced with the geodesic distance flSl , 
p4| . In | [25| , a method is proposed for clustering data lying 
on a manifold, which extends the graph-based semi-supervised 
learning algorithm in | [26| to a setting with unlabeled data. The 
diffusion matrix that diffuses known class labels to unlabeled 
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Fig. 1: PCA basis vectors computed with data sampled from a neighborhood on a manifold. In (a), the two most significant principal directions correspond to tangent directions 
and PCA computes a local model coherent with the manifold geometry. In (b), PCA fails to recover the tangent space as the manifold bends over itself and the neighborhood size 
is not selected properly. In (c), as the curvature component is stronger than the tangential components, the subspace spanned by the two most significant PCA basis vectors again 
fails to approximate the tangent space. 


data in | [26| is interpreted as a diffusion kernel in | [25| , 
which is then used for determining the similarity between data 
samples to obtain clusters. The works in p7| , also use 
the geodesic distance as a dissimilarity measure. They propose 
methods for embedding the manifold into the tangent spaces of 
some selected reference points and perform a fast approximate 
nearest neighbor search on the space of embedding. 

While the above algorithms consider all data samples to 
lie on a single manifold, several other methods model low¬ 
dimensional data as samples from multiple manifolds and 
study the determination of these manifolds. An expectation 
maximization approach is employed in | [29| to partition the 
data into manifolds. The points on each manifold are then 
embedded into a lower-dimensional domain. The method in 
p0| computes a sparse representation of each data sample in 
terms of other samples, where high coefficients are encouraged 
for nearby samples. Once the sparse coefficients are computed, 
data is grouped into manifolds simply with spectral clustering. 
The method in extends several popular nonlinear dimen¬ 
sionality reduction algorithms to the Riemannian setting by 
replacing the Euclidean distance with the Riemannian distance. 
It is then shown that, if most data connections lie within the 
manifolds rather than between them, the proposed Riemannian 
extensions yield clusters corresponding to different manifolds. 

Finally, the generation of overlapping clusters in GOC 
is also linked to soft clustering p2| . Rather than strictly 
partitioning the data into a set of disjoint groups, a membership 
score is computed between each data sample and each cluster 
center in soft clustering. The cluster centers are then updated 
by weighing the samples according to the membership scores. 
In p^ , a manifold extension of soft clustering is proposed, 
where the membership scores are computed with a geodesic 
kernel instead of the Euclidean distance. 

III. Rationale and Problem Formulation 

In patch-based image processing, one often would like 
to develop tools that can capture the common structures 
inherently present in patches and use this information for the 
efficient treatment of images. One important example is the 
invariance to geometric transformations. In practical image 
formation scenarios, different regions of the image are likely 
to observe the same structure, exposed, however, to different 


geometric transformations in different parts of the image 
plane. While most patch-based methods inherently achieve 
invariance to translations as they extract patches from the 
image over sliding windows, more complex transformations 
such as rotations and scale changes are more difficult to handle 
in evaluating the structural similarities between patches. In 
addition to geometric transformation models, structural simi¬ 
larities between image patches may be stemming from many 
other low-dimensional, possibly parametrizable patch models 
as well. In eg several parametrizable patch manifold models 
are explored such as oscillating textures and cartoon images. In 
the treatment or reconstruction of image patches, local models 
computed from patches sharing the same structure reflect the 
local geometry of the patch manifold, while the comparison 
of patch similarities based on Euclidean distance does not 
necessarily achieve this. In this paper, we propose similarity 
assessment strategies that better take structural similarities 
into account than the simple Euclidean distance in image 
reconstruction. 

Given observed measurements y, the ill-posed inverse prob¬ 
lem can be generally formulated in a Banach space as 

y = 0x + zy (1) 

where 0 is a bounded operator, x is an unknown data point 
and u is an error term. In image restoration y is the vectorized 
form of an observed image, 0 is a degradation matrix, x is 
the vectorized form of the original image, and u is an additive 
noise vector. There are infinitely many possible data points x 
that explain y; however, image restoration algorithms aim to 
reconstruct the original image x from the given measurements 
y, often by using some additional assumptions on x. 

In image restoration with sparse representations, x can be 
estimated by minimizing the cost function 

a = argmin |||y - 0$ oq ;||2 + A ||q:||i| (2) 

where ^ is a dictionary, a is the sparse representation of x in 
and A > 0 is a regularization parameter. It is common to 
reconstruct images patch by patch and model the patches of 
X as sparsely representable in Representing the extraction 
of the j-th patch Xj of x with a matrix multiplication as 
Xj = i?jX, the reconstruction of the overall image x can be 
represented via the operator o as shown in ||T|, Q. If the 
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dictionary ^ is well-chosen, one can efficiently model the 
data points x using their sparse representations in Once the 
sparse coefficient vector a is estimated, one can reconstruct 
the image x as 

X = ^ o a. (3) 

While a global model is considered in the above problem, 
several works such as Q. 0. p4| propose to reconstruct 
image patches based on sparse representations in local models. 
In this case, one aims to reconstruct the j-th patch Xj of 
the unknown image x from its degraded observation yj by 
selecting a local model that is suitable for yj . The problem in 
0 is then reformulated as 

aj = argmin hyj - Q^jajWl + A ||aj||i| (4) 

aj 

where yj is the j-th patch from the observed image y, is 
a local (PCA) basis chosen for the reconstruction of yj, and 
aj is the coefficient vector. The unknown patch Xj is then 
reconstructed as xj = ^j^j. The optimization problem in Q 
forces the coefficient vector aj to be sparse. Therefore, the 
accuracy of the reconstructed patch Xj in approximating the 
unknown patch Xj depends on the reliability of the local basis 
i.e., whether signals are indeed sparsely representable in 

The main idea proposed in this paper is to take into account 
the manifold structure underlying the data when choosing a 
neighborhood of training data points to learn a local basis. 
Our purpose is to develop a dissimilarity measure that is better 
suited to the local geometry of the data than the Euclidean dis¬ 
tance and also to make the neighborhood selection procedure 
as adaptive as possible to the test samples to be reconstructed. 

Let V = {di}^^ be a set of m training data points di G 
lying on a manifold A4 and let y = {yjj^i be a set of M test 
data points yj G As for the image reconstruction problem 
in (0, each test data point yj corresponds to a degraded image 
patch, and the training data points in V are used to learn the 
local bases ^j. The test samples yj are not expected to lie 
on the patch manifold A4 formed by the training samples; 
however, one can assume yj to be close to M. unless the 
image degradation is very severe. 

We then study the following problem. Given an observation 
yj G iy of an unknown image patch Xj, we would like to 
select a subset S CV of training samples such that the PCA 
basis ^j computed from S minimizes the reconstruction error 
\\xj — Xj\\, where the unknown patch Xj is reconstructed as 
Xj = ^jC^j, and the sparse coefficient vector is given by 

aj = argmin |||yj - Q^jajWl + A ||aj||d . (5) 

Since the nondeformed sample Xj is not known, it is clearly 
not possible to solve this problem directly. In this work, we 
propose some constructive solutions to guide the selection of 
S by assuming that yj lies close to Af. As the manifold M 
is not known analytically, we capture the manifold structure 
of training data V by building a similarity graph whose nodes 
and edges represent the data points and the affinities between 
them. In Sections [Iv] and lV| we describe the AGNN and the 
GOC methods, which respectively propose an adaptive and 



Fig. 2: Illustration of AGNN. The affinity between yj and di is ai, and the affinity 
between di and di is a*^. The intermediate node di contributes by the product aia*i 
to the overall affinity between yj and di. The sample is just another intermediate 
node like di. Summing the affinities via all possible intermediate nodes (i.e., all training 
samples), the overall affinity is obtained as in j^. 


a nonadaptive solution for training subset selection for local 
basis learning from the similarity graph. 


IV. Adaptive Geometry-Driven Nearest Neighbor 
Search 


In this section, we present the Adaptive Geometry-driven 
Nearest Neighbor search (AGNN) strategy for selecting the 
nearest neighbors of each test data point within the training 
data points with respect to an intrinsic manifold structure. 
Our subset selection method builds on the RGC algorithm 
P^ , which targets the clustering of data with respect to 
the underlying manifold. The RGC method seeks a globally 
consistent affinity matrix that is the same as its diffused version 
with respect to the underlying graph topology. However, the 
RGC method focuses only on the initially available training 
samples and does not provide a means of handling initially 
unavailable test samples. We thus present an out-of-sample 
generalization of RGC and propose a strategy to compute and 
diffuse the affinities between the test sample and all training 
samples in a way that is consistent with the data manifold. 

In the RGC algorithm, given a set of data points V, an 
affinity matrix A = (an) is first computed. The elements an 
of A measure the similarity between the data points di and di . 
A common similarity measure is the Gaussian kernel 


an = exp 


( nci^ j 


( 6 ) 


where H-H denotes the i 2 -^orm on and ci is a constant. 
Then, the initial affinities are updated with respect to the 
underlying manifold as follows. The affinities are diffused by 
looking for an A matrix such that each row Ai of A maximizes 

AJ = diTgmdix{v^ Av). (7) 

V 

Since the maximization problem on the right hand side of 0 
is solved by an eigenvector of A, the method seeks an affinity 
matrix such that the similarities between the data sample di 
and all the other samples in V (given by the row Ai) are 
proportional to the diffused version of the similarities in Ai 
over the whole manifold via the product AAj ; i.e., an affinity 
matrix is searched such that Aj oc AAj. The optimization 
problem in Q is solved with an iterative procedure based 
on a game theoretical approach to obtain a diffused affinity 
matrix A*. The diffusion of the affinities are constrained to 
the s nearest neighbors of each point di. 
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Algorithm 1 Adaptive Geometry-driven Nearest Neighbor 
search (AGNN) 

1: Input: 

T> = \ Set of training samples 

Vj E y: Test sample 

Cl, C 2 , K,: Algorithm parameters 

2: AGNN Algorithm: 

3: Form affinity matrix A of training samples with respect to 

4: Diffuse the affinities in A to obtain A* as proposed in the RGC method |18| . 

5: Initialize the affinity vector a between test sample yj and the training samples as 

in {8|. 

6: Diffuse the affinities in a to obtain a* with respect to jlO) . 

7: Determine set S of nearest neighbors of yj by selecting the training samples with 
the highest affinities as in (TT}. 

8: Output: 

S : Set of nearest neighbors of yj in T>. 


In our AGNN method, we first compute and diffuse the 
affinities of training samples in V as proposed in |T^ . This 
gives us a similarity measure coherent with the global geom¬ 
etry of the manifold. Meanwhile, unlike in RGC, our main 
purpose is to select a subset S CV of training samples for a 
given test sample yj e y. We thus need a tool for generalizing 
the above approach for test samples. 

We propose to compute the affinities between yj and V by 
employing A* as follows. Given a test data point yj G y, we 
first compute an initial affinity vector a whose i-th entry 


ai = exp 


I nci^ I 


( 8 ) 


measures the similarity between yj and the training sample 
di. We then update the affinity vector as follows. Denoting 
the entries of the diffused affinity matrix A* by a*^, first the 
product ahai should give the component of the overall affinity 
between yj and di that is obtained through the sample df. 
if there is a sample di that has a high affinity with both 
di and yj, this means that the affinity between di and yj 
should also be high due to the connection established via the 
intermediate node di (see the illustration in Figure [^. Note 
that the formulation in 0 also relies on the same idea. We 
thus update the affinity vector a such that its i-th entry 
becomes proportional to 

m 

1 = 1 

i.e., the total affinity between samples di and yj obtained 
through all nodes di in the training data graph. This suggests 
that the initial affinities in the vector a should be updated as 
A* a, which corresponds to the diffusion of the affinities on 
the graph. Repeating this diffusion process k, times, we get the 
diffused affinities of the test sample as 


= (A*)'"a 


( 10 ) 


where a* gives the final diffused affinity between yj and di. 
This generalizes the idea in Q to initially unavailable data 
samples; and hence, provides an out-of-sample extension of 
the diffusion approach in RGC. The parameter n. should be 
chosen in a way to permit a sufficient diffusion of the affinities. 
However, it should not be too large in order not to diverge too 
much from the initial affinities in a. In our experiments we 
have observed that n = 2 gives good results in general. 


Once the affinities a* are computed, the subset S consisting 
of the nearest neighbors of yj can be obtained as the samples 
in V whose affinities to yj are higher than a threshold 


S = {di e.T) \ al > C 2 maxa^} 


( 11 ) 


where 0 < C 2 < 1. The samples in S are then used for learning 
a PCA basis to reconstruct yj. The threshold C 2 should be 
chosen sufficiently high to select only the similar patches to 
the reference patch, however, it should not be selected too 
high in order to have sufficiently many neighbors necessary 
for computing a basis. If S contains too few samples, the 
threshold C 2 can be adapted to increase the number of samples 
or a sufficient number of points with highest affinities can 
be directly included in S. The proposed AGNN method for 
determining training subsets gets around the problem depicted 
in Figure |l(b) since points lying at different sides of a 
manifold twisting onto itself have a small diffused affinity 
and are not included in the same subset. A summary of the 
proposed AGNN method is given in Algorithm 


V. Geometry-Driven Overlapping Clusters 


As we will see in Section VI the AGNN method presented 
in Section |lv| is efficient in terms of image reconstruction 
performance. However, it may have a high computational 
complexity and considerable memory requirements in settings 
with a large training set V, as the size of the affinity matrix 
grows quadratically with the number of training samples and 
the subset selection is adaptive (repeated for each test sample). 
For this reason, we propose in this section the Geometry- 
driven Overlapping Clusters (GOC) method, which provides 
a computationally less complex solution for obtaining the 
nearest neighbors of test samples. 

The GOC algorithm computes a collection {Sk}^^i of 
subsets Sk dV of the training data set, which are to be used 
in local basis computation. Contrary to the AGNN method, 
the subsets Sk C V are determined only using the training 
data and are not adapted to the test samples. However, the 
number C of subsets should then be sufficiently large to have 
the desired adaptivity for capturing arbitrary local variations. 
Due to the large number of subsets, Sk are not disjoint in 
general; hence, can be regarded as overlapping clusters. In 
the following, we first describe our method for forming the 
clusters and then propose a strategy to select some parameters 
that determine the size and the structure of the clusters. 

Given the number of clusters C to be formed, we first 
determine the central data point G D of each cluster Sk- 
In our implementation, we achieve this by first clustering V 
with the K-means algorithm, and then choosing each /i/c as 
the point in V that has the smallest Euclidean distance to the 
center of the k-th cluster given by K-means. 

The training data points yk are used as the kickoff for the 
formation of the clusters Sk- Given the central sample /i/^, 
the cluster Sk is formed iteratively with the GOC algorithm 
illustrated in Figure as follows. We first initialize Sk as 

Sk = 


( 12 ) 
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Fig. 3: Illustration of the GOC algorithm. The cluster Sk around the central sample fik 
is formed gradually. Sk is initialized with containing the K nearest neighbors of 
Hk = 3 in the illustration). Then in each iteration I, is expanded by adding the 
nearest neighbors of recently added samples. 


where MK{l^k) denotes the set of the if-nearest neighbors of 
Ilk in V with respect to the Euclidean distance. Then in each 
iteration /, we update the cluster S\. as 

Si = Sl-^yj U MM (13) 

dies^-^ 

by including all samples in the previous iteration as well as 
their if-nearest neighbors. Hence, the clusters are gradually 
expanded by following the nearest neighborhood connections 
on the data graph. This procedure is repeated for L iterations 
so that the final set of clusters is given by 

{Sk}k=l = {^k}k=l' (14) 

The expansion of the clusters is in a similar spirit to the affinity 
diffusion principle of AGNN; however, is computationally 
much less complex. 

In the simple strategy presented in this section, we have 
two important parameters to set, which essentially infiuence 
the performance of learning: the number of iterations L and 
the number of samples K in each small neighborhood. In the 
following, we propose an algorithm to adaptively set these 
parameters based on the local geometry of data. Our method 
is based on the observation that the samples in each cluster 
will eventually be used to learn a local subspace that provides 
an approximation of the local tangent space of the manifold. 
Therefore, Sk should lie close to a low-dimensional subspace 
in so that nearby test samples can be assumed to have a 
sparse representation in the basis ^k computed from Sk- We 
characterize the concentration of the samples in Sk around a 
low-dimensional subspace by the decay of the coefficients of 
the samples in the local PC A basis. 

We omit the cluster index k for a moment to simplify 
the notation and consider the formation of a certain cluster 
S = Sk- With a slight abuse of notation, let 5'^’^ stand for the 
cluster S that is computed by the algorithm described above 
with parameters L and K. Let ^ = [0i ... 0^] be the PCA 
basis computed with the samples in S, where the principal 
vectors 01 ,...,^^ ^ are sorted with respect to the 
decreasing order of the absolute values of their corresponding 
eigenvalues. For a training point di G S, let di = di — rjs 
denote the shifted version of di where r]s = 


Algorithm 2 Geometry-driven Overlapping Clusters (GOC) 

1: Input: 

^ = {di}^-^'. Set of training samples 
C\ Number of clusters 
C 3 : Algorithm parameter 

2: GOC Algorithm: 

3: Determine cluster centers /Xfc of all C clusters (possibly with the K-means algo¬ 
rithm). 

4: for /c = 1, • • • , C do 

5: Fix parameter L' = Lq at an initial value Lq. 

6: for K' = 1, ■ ■ ■ , Kmax do 

7: Form cluster Sk = as described in 

8: Evaluate decay rate function I{Lo, K') given in jib) . 

9: end for 

10: Set K as the K' value that minimizes i(Lo, K'). 

11. for L =!,••• ^ L max do 

12: Form cluster Sk = S^ as described in fI^- |T4l . 

13: Evaluate decay rate function I{L', K) given by (16) . 

14: end for 

15: Set L as the L' value that mi ni mi zes I{L', K). 

16: Determine cluster Sk as S^’^ with the optimized parameters. 

17: end for 
18: Output: 

: Set of overlapping clusters in D. 


is the centroid of cluster S. We define 

L 

I{L,K) = min{^ | E E 

9=1 

n 0^1 

<?=i dies^^^ 

which gives the smallest number of principal vectors to 
generate a subspace that captures a given proportion C 3 of 
the total energy of the samples in S, where 0 < C 3 < 1 . We 
propose to set the parameters L, K by minimizing the function 
which gives a measure of the concentration of the 
energy of S around a low-dimensional subspace. However, 
in the case that S contains m < n samples where n is the 
dimension of the ambient space, the subspace spanned by 
the first m — 1 principal vectors always captures all of the 
energy in S\ therefore I{L^K) takes a relatively small value; 
i.e., I{L,K) < m — 1. In order not to bias the algorithm 
towards reducing the size of the clusters as a result of this, a 
normalization of the function I{L,K) is required. We define 

I{L,K) = I{L,K)/mm{ \S^’^\ - 1, n} (16) 

where | • | denotes the cardinality of a set. The denomina¬ 
tor min{ 1 5'^’^ I — 1 , n} of the above expression gives the 
maximum possible value of /(L, AT) in cluster 5'^’^. Hence, 
the normalization of the coefficient decay function by its 
maximum value prevents the bias towards small clusters. 

We can finally formulate the selection of L, iT as (L, if) = 
arg min(x,',K')GA where A is a bounded parameter 

domain. This optimization problem is not easy to solve exactly. 
One can possibly evaluate the values of I{L^K) on a two- 
dimensional grid in the parameter domain. However, in order 
to reduce the computation cost, we approximately minimize 
the objective by optimizing one of the parameters and fixing 
the other in each iteration. We first fix the number of iterations 
L at an initial value and optimize the number of neighbors K. 
Then, updating and fixing if, we optimize L. 

The computation of the parameters L and K with the above 
procedure determines the clusters as in ( p^ . The samples in 
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each cluster Sk are then used for computing a local basis 
The proposed GOC method is summarized in Algorithm 
Since the proposed GOC method determines the clusters not 
only with respect to the connectivity of the data samples on 
the graph, but also by adjusting the size of the clusters with 
respect to the local geometry, it provides a solution for both 
of the problems described in Figures |l(b)| and 

In the proposed GOC method, contrary to AGNN, we need 
to define a strategy to select the PCA basis that best fits a given 
test patch. Given a test patch i/j , we propose to select a basis 
by taking into account the distance between yj and the 
centroid yk of the cluster Sk (corresponding to T>/c), as well 
as the agreement between yj and the principal directions in ^k- 
Let = [ 01 ... 0r] denote the submatrix of ^k consisting 
of the first r principal vectors, which give the directions that 
determine the main orientation of the cluster. We then choose 
the basis ^k that minimizes 


1(c) 
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Fig. 4: Two of the reference patches and their rotated versions used in the experiment 



Fig. 5: Percentage of patches correctly included in the clusters 


k = arg min 
k' 
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(17) 

where 7 > 0 is a weight parameter. While the first term above 
minimizes the distance to the centroid of the cluster, the second 
term maximizes the correlation between the relative patch 
position yj — yk' and the most significant principal directions. 
Once the basis index k is determined as above, the test patch 
yj is reconstructed based on a sparse representation in ^k- 


VI. Experiments 

We verify the performance of our proposed methods with 
extensive experiments on image restoration based on sparse 
representations. In Section [VI-A| we first present an experiment 
where we evaluate the performance of the proposed neighbor¬ 
hood selection strategies in capturing the structural similarities 
of images. Then in Sections |VI-B[ |VI-C[ and |VI-D| we test 
our algorithms respectively in superresolution, deblurring, and 
denoising applications. 


A. Transformation-invariant patch similarity analysis 

Natural images often contain different observations of the 
same structure in different regions of the image. Patches that 
share a common structure may be generated from the same 
reference pattern with respect to a transformation model that 
can possibly be parameterized with a few parameters. One 
example to parametrizable transformation models is geometric 
transformations. In this section, we evaluate the performance 
of the proposed AGNN strategy in capturing structural similar¬ 
ities between image patches in a transformation-invariant way. 
We generate a collection of patches of size 10 x 10 pixels, by 
taking a small set of reference patches and applying geometric 
transformations consisting of a rotation with different angles 
to each reference patch to obtain a set of geometrically 
transformed versions of it. Figure shows two reference 
patches and some of their rotated versions. The data set used in 
the experiment is generated from 10 reference patches, which 
are rotated at intervals of 5 degrees. 


In order to evaluate the performance of transformation- 
invariant similarity assessment, we look for the nearest neigh¬ 
bors of each patch in the whole collection and identify the 
“correct” neighbors as the ones sharing the same structure, 
i.e., the patches generated from the same reference patch. 
Three nearest neighbor selection strategies are tested in the 
experiment, which are AGNN, neighbor selection with respect 
to Euclidean distance, and K-means clustering. In AGNN, the 
neighborhood size that gives the best algorithm performance 
is used. The Euclidean distance uses the same neighborhood 
size as AGNN, and the number of clusters in K-means is set 
as the true number of clusters, i.e., the number of reference 
patches generating the data set. The correct clustering rates are 
shown in Figure which are the percentage of patches that are 
correctly present in a cluster (each neighborhood is considered 
as a cluster in AGNN and Euclidean distance). The horizontal 
axis shows the number of clusters (i.e., number of reference 
patches) used in different repetitions of the experiment. It 
can be observed that the AGNN method yields the best 
transformation-invariant similarity assessment performance. 
Contrary to methods based on simple Euclidean distance, 
AGNN measures the similarity of two patches by tracing all 
paths on the manifold joining them. Therefore, it is capable 
of following the gradual transformations of structures on the 
patch manifold and thus identifying structural similarities of 
patches in a transformation-invariant manner. 


B. Image super re solution 

In this section, we demonstrate the benefits of our neighbor¬ 
hood selection strategies in the context of the NCSR algorithm 
Q, which leads to state-of-the-art performance in image 
superresolution. 

The NCSR algorithm Q is an image restoration method 
that reconstructs image patches by selecting a model among 
a set of local PCA bases. This strategy exploits the image 
nonlocal self-similarity to obtain estimates of the sparse coding 
coefficients of the observed image. The method first clusters 
training patches with the K-means algorithm and then adopts 
the adaptive sparse domain selection strategy proposed in 
























Fig. 6: Test images for super-resolution: Butterfly, Bike, Hat, Plants, Girl, Parrot, 
Parthenon, Raccoon, Leaves, Flower. 


to learn a local PCA basis for each cluster from the 
estimated high-resolution (HR) images. After the patches are 
coded, the NCSR objective function is optimized with the 
Iterative Shrinkage Thresholding (1ST) algorithm proposed in 
| [35| . The clustering of training patches with the K-means 
algorithm in Q is based on adopting the Euclidean distance 
as a dissimilarity measure. The purpose of our experiments 
is then to show that the proposed geometry-based nearest 
neighbor selection methods can be used for improving the 
performance of an image reconstruction algorithm such as 
NCSR. 

We now describe the details of our experimental setting 
for the superresolution problem. In the inverse problem y = 
0 x + z/ in Q, X and y denote respectively the lexicographical 
representations of the unknown image X and the degraded 
image Y. The degradation matrix 0 = DH is composed of a 
down-sampling operator D with a scale factor ofq = 3 and a 
Gaussian filter H of size 7x7 with a standard deviation of 1.6, 
and u is an additive noise. We aim to recover the unknown 
image vector x from the observed image vector y. We evaluate 
the proposed algorithms on the 10 images presented in Figure 

which differ in their frequency characteristics and content. 
For color images, we apply the single image SR algorithm 
only on the luminance channel and we compute the PSNR 
and SSIM | [36| only on the luminance channel for coherence. 
Besides PSNR and SSIM, the visual quality of the images is 
also used as a comparison metric. 

Overlapping patches of size 6 x 6 are used in the exper¬ 
iments. The original NCSR algorithm initializes the training 
set V by extracting patches from several images in the scale 
space of the HR image. However, in our implementation we 
initialize the set of training patches by extracting them only 
from the low-resolution image; i.e., the m initial training 
patches di G in V = are extracted from the 

observed low-resolution (LR) image vector y. We learn online 
PCA bases using the training patches in V with the proposed 
AGNN and GOC methods. In the original NCSR method, in 
every P iterations of the 1ST algorithm, the training set V 
is updated by extracting the training patches from the current 
version of the reconstructed image x and the PCA bases are 
updated as well by repeating the neighborhood selection with 
the updated training data. In our experiments, we use the same 
training patches V for the whole algorithm. 

In Section | VI-B1 1 we evaluate our methods AGNN and 
GOC by comparing their performance to some other clustering 


or nearest neighbor selection strategies in superresolution. 
In Section |VI-B2[ we provide comparative experiments with 
several widely used superresolution algorithms and show that 
our proposed manifold-based neighborhood selection tech¬ 
niques can be used for improving the state of the art in 
superresolution. 

In Sections | VI-B 1 and VI-B2 the results for the algorithms 
in comparison are obtained by using the software packages 
made publicly available by the corresponding authors Q 

1) Performance Evaluation of AGNN and GOC: We com¬ 
pare the proposed AGNN and GOC methods with 4 dif¬ 
ferent clustering algorithms; namely, the K-means algorithm 
(Kmeans), Fuzzy C-means clustering algorithm (FCM) p2| . 
Spectral Clustering (SC) p0| . Replicator Graph Clustering 
(RGC) fT^ ; and also with K-NN search using geodesic 
distance (GeoD). Among the clustering methods, Kmeans 
and FCM employ the Euclidean distance as a dissimilarity 
measure, while SC and RGC are graph-based methods that 
consider the manifold structure of data. When testing these 
four methods, we cluster the training patches and compute a 
PCA basis for each cluster. Then, given a test patch, the basis 
of the cluster whose centroid has the smallest distance to the 
test patch is selected as done in the original NCSR algorithm 
where K-means is used. In the GeoD method, each test patch is 
reconstructed with the PCA basis computed from its nearest 
neighbors with respect to the geodesic distance numerically 
computed with Dijkstra’s algorithm p7| . The idea of nearest 
neighbor selection with respect to the geodesic distance is also 
in the core of the methods proposed in p7| and 1^ . Note 
that the four reference clustering methods and GOC provide 
nonadaptive solutions for training subset selection, while the 
GeoD and the AGNN methods are adaptive. 

The parameters of the AGNN algorithm are set as s = 35 
(number of nearest neighbors in the diffusion stage of RGC 
GD) , n = 2 (number of iterations for diffusing the affinity 
matrix), ci = 10 (Gaussian kernel scale), and C 2 = 0.9 
(affinity threshold). The parameters of the GOC algorithm 
are set as C = 64 (number of clusters), cs = 0.5 (threshold 
defining the decay function), 7 = 150, and r = 8 (parameters 
for selecting a PCA basis for each test patch). The number of 
clusters in the other four clustering methods in comparison are 
also set to the same value as C = 64. The size of the clusters 
with the FCM algorithm are selected roughly the same as the 
cluster sizes computed with K-means. The total number of 
iterations and the number of PCA basis updates are chosen as 
1000 and 4 in the NCSR algorithm. All the general parameters 
for the NCSR algorithm are selected as Dong et al. Q. In 
this way, we can maintain consistency in the comparison of 
the methods related to NCSR algorithm. 

We evaluate the GOC algorithm in three different settings. 
In the first setting the cluster size parameters L and K 
are estimated adaptively for each cluster with the strategy 
proposed in Algorithm which is denoted as aGOC. In the 
second setting, denoted avGOC, the parameters L and K are 
not adapted to each cluster; all clusters are formed with the 


^We would like to thank the authors of (^, ( 2 ), (^, (T^ , |2ol , and 
El for making their software packages publicly available. 
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(a) LR image (b) Original HR image (c) NCSR-Kmeans (d) NCSR-AGNN (e) NCSR-GOC 



(f) Original HR close-up (g) NCSR-Kmeans close-up (h) NCSR-AGNN close-up (i) NCSR-GOC close-up 


Fig. 7: Comparison of SR results (x3). It can be observed that NCSR-AGNN and NCSR-GOC reconstruct edges with a higher contrast than NCSR-Kmeans. Artifacts visible with 
NCSR-Kmeans (e.g., the oscillatory phantom bands perpendicular to the black stripes on the butterfly’s wing) are significantly reduced with NCSR-AGNN and NCSR-GOC. 



(a) LR image (b) Original HR image (c) NCSR-Kmeans (d) NCSR-AGNN (e) NCSR-GOC 



(f) Original HR close-up (g) NCSR-Kmeans close-up (h) NCSR-AGNN close-up (i) NCSR-GOC close-up 


Fig. 8: Comparison of SR results (x3). NCSR-Kmeans produces artifacts such as the checkerboard-like noise patterns visible on plain regions of the cap, which are prevented by 
NCSR-AGNN or NCSR-GOC. 


same parameter values, where L and K are computed by 
minimizing the average value of coefficient decay function 
I{L,K) over all clusters of the same image. The parameters 
are thus adapted to the images, but not to the individual 
clusters of patches of an image. Finally, in the third setting, 
denoted mGOC, the parameters L and K are manually entered 
and used for all clusters of the same image. The parameter 
values provided to the algorithm for each image are set as 
the best values obtained with an exhaustive search. Therefore, 
mGOC can be considered as an oracle setting. 

The results are presented in Figure Figure and Table 
|l| Figures [7] and provide a visual comparison between 
the image reconstruction qualities obtained with the K-means 
clustering algorithm and the proposed AGNN and GOC meth¬ 
ods for the Butterfly and the Hat images. It is observed 
that AGNN and GOC produce sharper edges than K-means. 


Moreover, the visual artifacts produced by K-means such as 
the phantom perpendicular bands on the black stripes of the 
butterfly and the checkerboard-like noise patterns on the cap 
are significantly reduced with AGNN and GOC. The efficiency 
of the proposed methods for removing these artifacts can 
be explained as follows. When image patches are clustered 
with the K-means algorithm, the similarity between patches 
is measured with the Euclidean distance. Therefore, when 
reconstructing a test patch, the algorithm tends to use a basis 
computed with patches that have similar intensity values. The 
nonuniformity of the pixel intensities along the black stripes 
of the LR Butterfiy image thus propagates to the reconstructed 
HR image as well, which produces the phantom bands on the 
wing (due to the too low resolution, the black stripes on the LR 
image contain periodically appearing clear pixels contaminated 
by the yellow plain regions on the wing). Similarly, in the Hat 
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TABLE I: PSNR (top row, in dB) and SSIM (bottom row) results for the luminance components of super-resolved HR images for different clustering or neighborhood selection 
approaches: Spectral Clustering (SC) Fuzzy C-means clustering algorithm (FCM) K-means clustering (Kmeans); Replicator Graph Clustering (RGC) f^; kNN search 


Images 

Butterfly 

Bike 

Hat 

Plants 

Girl 

Parrot 

Parthenon 

Raccoon 

Leaves 

Flower 

Average 

SC 

28.15 

24.73 

31.28 

33.98 

33.65 

30.45 

27.19 

29.24 

27.50 

29.45 

29.56 

0.9193 

0.8026 

0.8723 

0.9198 

0.8255 

0.9170 

0.7509 

0.7659 

0.9242 

0.8567 

0.8554 

FCM 

28.20 

24.76 

31.25 

33.99 

33.65 

30.47 

27.25 

29.25 

27.68 

29.50 

29.60 

0.9205 

0.8040 

0.8726 

0.9205 

0.8256 

0.9174 

0.7531 

0.7663 

0.9271 

0.8575 

0.8565 

Kmeans 

28.14 

24.79 

31.31 

34.07 

33.64 

30.53 

27.20 

29.28 

27.67 

29.47 

29.61 

0.9204 

0.8050 

0.8730 

0.9213 

0.8254 

0.9178 

0.7517 

0.7668 

0.9265 

0.8567 

0.8565 

RGC [Is] 

28.45 

24.80 

31.37 

34.20 

33.65 

30.57 

27.22 

29.27 

27.90 

29.50 

29.69 

0.9234 

0.8061 

0.8739 

0.9219 

0.8254 

0.9181 

0.7525 

0.7658 

0.9317 

0.8576 

0.8576 

GeoD 

28.61 

24.82 

31.42 

34.16 

33.63 

30.44 

27.24 

29.25 

27.98 

29.54 

29.71 

0.9257 

0.8070 

0.8746 

0.9219 

0.8250 

0.9178 

0.7530 

0.7650 

0.9323 

0.8587 

0.8581 

avGOC 

28.34 

24.85 

31.42 

34.17 

33.66 

30.68 

27.23 

29.28 

27.89 

29.55 

29.71 

0.9222 

0.8076 

0.8747 

0.9224 

0.8258 

0.9191 

0.7528 

0.7668 

0.9317 

0.8591 

0.8582 

aGOC 

28.46 

24.85 

31.44 

34.18 

33.65 

30.63 

27.23 

29.27 

27.92 

29.54 

29.72 

0.9239 

0.8082 

0.8744 

0.9227 

0.8257 

0.9187 

0.7530 

0.7663 

0.9324 

0.8588 

0.8584 

mGOC 

28.54 

24.90 

31.43 

34.20 

33.67 

30.71 

27.25 

29.28 

27.95 

29.55 

29.75 

0.9251 

0.8085 

0.8748 

0.9222 

0.8261 

0.9192 

0.7530 

0.7671 

0.9324 

0.8593 

0.8588 

AGNN 

28.78 

24.87 

31.46 

34.16 

33.67 

30.60 

27.29 

29.26 

28.01 

29.61 

29.77 

0.9266 

0.8081 

0.8749 

0.9218 

0.8260 

0.9188 

0.7540 

0.7661 

0.9324 

0.8601 

0.8589 


image, the clusters used in learning a basis for reconstructing 
the plain regions on the cap contain also patches extracted 
from the wall, which have a similar intensity with the cap. 
This reproduces the shadowy patterns of the wall also on the 
cap. On the other hand, the AGNN method groups together 
patches that have a connection on the data graph. As the 
patches are extracted with overlapping windows shifting by 
one pixel, AGNN and GOC may have a stronger tendency than 
K-means for favoring patches from nearby or similar regions 
on the image that all share a common structure, which is also 
confirmed by the experiment in Section VI-A The proposed 
methods yield local bases better fitted to the characteristics of 
patches, therefore, less artifacts are observed. 


In Table |T| the performance of the compared clustering 
methods are measured with the PSNR and the SSIM metrics. 
Graph-based methods are generally seen to yield a better 
performance than methods based on Euclidean distance. This 
confirms the intuition that motivates our study; when selecting 
neighborhoods for learning local models, the geometry of the 
data should be respected. As far as the average performance is 
concerned, the AGNN method gives the highest reconstruction 
quality and is followed by the GOC method. The performance 
difference between AGNN and GOC can be justified with 
the fact that the training subset selection is adaptive to the 
test patches in AGNN, while GOC is a nonadaptive method 
that offers a less complex solution. In particular, with a non- 
optimized implementation of our algorithms, we have observed 
that GOC has roughly the same computation time as K-means, 
while the computation time of AGNN is around three times 
K-means and GOC in the tested images on an Intel Core i5 
2.6GHz under the Matlab R2015a programming environment, 
as shown in Table |nl] K-means and GOC in the tested images. 
After the proposed AGNN and GOC methods, GeoD gives the 


best average performance. While this adaptive method ensures 
a good reconstruction quality, it requires the computation of 
the geodesic distance between each test patch and all training 
patches. Therefore, it is computationally very complex. Al¬ 
though several works such as p7| and p8| provide solutions 
for fast approximations of the geodesic distance, we observe 
that in terms of reconstruction quality AGNN performs better 
than GeoD in most images. This suggests that using a globally 
consistent affinity measure optimized with respect to the entire 
graph topology provides a more refined and precise similarity 
metric than the geodesic distance, which only takes into 
account the shortest paths between samples. 

Concerning the performances of the clustering methods 
on the individual images, an important conclusion is that 
geometry-based methods yield a better performance especially 
for images that contain patches of rich texture. The AGNN 
and GOC methods provide a performance gain of respectively 
0.64 dB and 0.4 dB over K-means (used in the original NCSR 
method) for the Butterfiy image. Meanwhile, all clustering 
methods give similar reconstruction qualities for the Girl 
image. This discrepancy can be explained with the difference 
in the characteristics of the patch manifolds of these two 
images. The patches of the Butterfiy image contain high- 
frequency textures; therefore, the patch manifold has a large 
curvature (see, e.g., | [3^ for a study of the relation between 
the manifold curvature and the image characteristics). Conse¬ 
quently, the proposed methods adapted to the local geometry 
of the manifold perform better on this image. On the other 
hand, the Girl image mostly contains weakly textured low- 
frequency patches, which generate a rather fiat patch manifold 
of small curvature. The Euclidean distance is more reliable as 
a dissimilarity measure on fiat manifolds compared to curved 
manifolds as it gets closer to the geodesic distance. Hence, the 
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performance gain of geometry-based methods over K-means 
is much smaller on the Girl image compared to Butterfly. 

Next, the comparison of the three modes of the GOC 
algorithm shows that aGOC and avGOC yield reconstruction 
qualities that are close to that of the oracle method mGOC. 
This suggests that setting the parameters L and K with respect 
to the PC A coefficient decay rates as proposed in Algorithm]^ 
provides an efficient strategy for the automatic determination 
of cluster sizes. While the average performances of aGOC 
and avGOC are quite close, interestingly, aGOC performs 
better than avGOC on Butterfly and Leaves. Both of these 
two images contain patches of quite varying characteristics, 
e.g., highly textured regions formed by repetitive edges as 
well as weakly textured regions. As the structures of the 
patches change signiflcantly among different clusters in these 
images, optimizing the cluster size parameters individually for 
each cluster in aGOC has an advantage over using common 
parameters in avGOC. 

2) Improvements over the State of the Art in Superreso¬ 
lution: In this section, we present an experimental compar¬ 
ison of several popular superresolution algorithms; namely, 
the bicubic interpolation algorithm, ASDS 0 SPSR (^, 
and NCSR 0. We evaluate the performance of the NCSR 
algorithm under three different settings where the local bases 
are computed with K-means, AGNN, and GOC. The GOC 
method is used as in Algorithm (denoted as aGOC in the 
previous experiments). 

The experiments are conducted on the same images as in 
the previous set of experiments. The total number of iterations 
and the number of PCA basis updates of NCSR are selected 
respectively as 960 and 6 , while the other parameters are 
chosen as before. The results presented in Table |I^ show 
that the state of the art in superresolution is led by the 
NCSR method Q. The performance of NCSR is improved 
when it is coupled with the AGNN and GOC strategies for 
selecting local models. In Table the images are divided 
into two categories as those with high-frequency and low- 
frequency content. The average PSNR and SSIM metrics are 
reported in both groups. It can be observed that the advantage 
of the proposed neighborhood selection strategies over K- 
means is especially signiflcant for high-frequency images. 
In images with low-frequency content, K-means gives the 
same performance as the proposed methods. As the patch 
manifold gets flatter, clusters obtained with K-means and the 
proposed methods get similar. Hence, we may conclude that 
the proposed geometry-based neighborhood selection methods 
can be successfully used for improving the state of the art in 
image superresolution, whose efficacy is especially observable 
for sharp images rich in high-frequency texture. 


C. Image deblurring 

We now evaluate our method in the image deblurring 
application. Unlike the case in superresolution, the images 
to be deblurred have a normal resolution, which leads to a 
large number of patches for large images. In this case GOC 
has an advantage over AGNN in terms of complexity and 
memory requirements. Thus, it is more interesting to study 



Fig. 9: Test images for deblurring: Butterfly, Boats, Cameraman, House, Parrot, Lena, 
Barbara, Starfish, Peppers, Leaves. 


the performance of the GOC algorithm in deblurring. We 
compare GOC with the K-means clustering algorithm within 
the framework of the NCSR method jf^. The algorithms are 
tested on the images shown in Figure Two blurring kernels 
are used, which are a uniform blur kernel of size 9x9 pixels 
and a Gaussian blur kernel of standard deviation 1.6 pixels. 
Along with the blurring, the images are also corrupted with an 
additive white Gaussian noise of standard deviation s/2. The 
parameters of GOC are set as C = 64 (number of clusters), 
C 3 = 0.5 (threshold deflning the decay function), 7 = 150, 
and r = 8 (parameters for selecting a PCA basis for each test 
patch). All the general parameters for the NCSR algorithm are 
selected as Dong et al. |[T| in order to maintain the consistency. 

The PSNR and FSIM |[39l measures of the reconstruction 


qualities are presented in Table IV The results obtained with 


the image restoration algorithms FISTA (Portilla et al.) |40 
/q-SPAR (Irani et al.) |[4T|, IDD-BM3D (Danielyan et al.) |42 


and ASDS (Dong et alOQ reported in Q for the same experi¬ 
ments are also given for the purpose of comparison. The results 
show that the proposed GOC algorithm can be effectively used 
for improving the image reconstruction quality of the NCSR 
method in deblurring applications. The GOC method either 
outperforms the K-means clustering algorithm or yields a quite 
close performance when coupled with NCSR. Moreover, one 
can observe that the best average PSNR value is given by the 
proposed method, whose benefits are especially observable for 
images with signiflcant high-frequency components such as 
Butterfly, Cameraman, and Leaves. 


D. Image denoising 

We now evaluate our method in the image denoising appli¬ 
cation. Since the deformation of the patch manifold geometry 
due to noise poses a challenge on geometry-based similarity 
assessment between patches, we use the AGNN method in the 
experiments in this section, which usually has a better recon¬ 
struction quality than GOC. We compare AGNN with K-means 
within the framework of the NCSR method j[B. The algorithms 
are tested on the images shown in Figure The images are 
corrupted with additive white Gaussian noise at different noise 
levels with standard deviation a = [5 10 15 20 50 100]. 
The parameters of AGNN are set as s = 35 (number of nearest 
neighbors in the diffusion stage of RGC |T^), n = 2 (number 
of iterations for diffusing the affinities), ci = 10 (Gaussian 
kernel scale), and C 2 = 0.9 (affinity threshold). All the general 
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TABLE II: PSNR (top row, in dB) and SSIM (bottom row) results for the luminance components of super-resolved HR images for different superresolution algorithms: Bicubic 
Interpolation; SPSR (Peleg et al.) ASDS (Dong et al.) [^; NCSR (Dong et al.) NCSR with proposed GOC; NCSR with proposed AGNN. The methods are ordered 
according to the average PSNR values (from the lowest to the highest). 


Images 

Butterfly 

Bike 

Hat 

Plants 

Leaves 

Average 

Parrot 

Parthenon 

Raccoon 

Girl 

Elower 

Average 

Bicubic 

22.41 

21.77 

28.22 

29.69 

21.73 

24.76 

26.54 

25.20 

27.54 

31.65 

26.16 

27.42 

0.7705 

0.6299 

0.8056 

0.8286 

0.7302 

0.7530 

0.8493 

0.6528 

0.6737 

0.7671 

0.7295 

0.7345 

SPSR [ig 

26.74 

24.31 

30.84 

32.83 

25.84 

28.11 

29.68 

26.77 

29.00 

33.40 

28.89 

29.55 

0.8973 

0.7830 

0.8674 

0.9036 

0.8892 

0.8681 

0.9089 

0.7310 

0.7562 

0.8211 

0.8415 

0.8117 

ASDS (3 

27.34 

24.62 

30.93 

33.47 

26.80 

28.63 

30.00 

26.83 

29.24 

33.53 

29.19 

29.76 

0.9047 

0.7962 

0.8706 

0.9095 

0.9058 

0.8774 

0.9093 

0.7349 

0.7677 

0.8242 

0.8480 

0.8168 

NCSR [T] 

28.07 

24.74 

31.29 

34.05 

27.46 

29.12 

30.49 

27.18 

29.27 

33.66 

29.50 

30.02 

0.9156 

0.8031 

0.8704 

0.9188 

0.9219 

0.8860 

0.9147 

0.7510 

0.7707 

0.8276 

0.8563 

0.8241 

NCSR-GOC 

28.47 

24.85 

31.44 

34.16 

28.05 

29.39 

30.71 

27.23 

29.28 

33.65 

29.58 

30.09 

0.9241 

0.8084 

0.8747 

0.9232 

0.9339 

0.8929 

0.9192 

0.7526 

0.7666 

0.8257 

0.8600 

0.8248 

NCSR-AGNN 

28.81 

24.86 

31.47 

34.19 

28.06 

29.48 

30.60 

27.30 

29.27 

33.67 

29.60 

30.09 

0.9273 

0.8080 

0.8755 

0.9223 

0.9332 

0.8933 

0.9189 

0.7546 

0.7662 

0.8261 

0.8601 

0.8252 


TABLE III: Running times for the luminance components of super-resolved HR images for different super-resolution algorithms: NCSR (Dong et al.) NCSR with proposed 
GOC; NCSR with proposed AGNN. 


Images 

Butterfly 

Bike 

Hat 

Plants 

Leaves 

Parrot 

Parthenon 

Raccoon 

Girl 

Elower 

Average 

NCSR [T] 

261 

229 

213 

229 

233 

220 

481 

362 

213 

226 

267 

NCSR-GOC 

271 

266 

253 

261 

278 

256 

518 

383 

246 

264 

299 

NCSR-AGNN 

960 

1039 

467 

578 

1146 

505 

2541 

1637 

416 

830 

1012 


TABLE IV: PSNR (top row, in dB) and ESIM (bottom row) results for the luminance components of deblurred images for different debl urring algorithms for un ifor m blur kernel 
and Gaussian blur kernel of standard deviation 1.6 pixels: NCSR (Dong et al.) NCSR with proposed GOC; EISTA (Portilla et al.) j^; Zq-SPAR (Irani et al.) j^; IDD-BM3D 


Images 

Butterfly 

Boats 

C. Man 

House 

Parrot 

Lena 

Barbara 

Starfish 

Peppers 

Leaves 

Average 


Uniform 

EISTA 

28.37 

29.04 

26.82 

31.99 

29.11 

28.33 

25.75 

27.75 

28.43 

26.49 

28.21 

0.9119 

0.8858 

0.8627 

0.9017 

0.9002 

0.8798 

0.8375 

0.8775 

0.8813 

0.8958 

0.8834 


27.10 

29.86 

26.97 

32.98 

29.34 

28.72 

26.42 

28.11 

28.66 

26.30 

28.44 

Zo-SPAR (4T] 

0.8879 

0.9094 

0.8689 

0.9225 

0.9262 

0.9063 

0.8691 

0.8951 

0.9066 

0.8776 

0.8970 

ASDS (3 

28.70 

30.80 

28.08 

34.03 

31.22 

29.92 

27.86 

29.72 

29.48 

28.59 

29.84 

0.9053 

0.9236 

0.8950 

0.9337 

0.9306 

0.9256 

0.9088 

0.9208 

0.9203 

0.9075 

0.9171 

IDD-BM3D 

29.21 

31.20 

28.56 

34.44 

31.06 

29.70 

27.98 

29.48 

29.62 

29.38 

30.06 

0.9287 

0.9304 

0.9007 

0.9369 

0.9364 

0.9197 

0.9014 

0.9167 

0.9200 

0.9295 

0.9220 

NCSR [T] 

29.73 

31.04 

28.61 

34.26 

31.98 

29.95 

28.07 

30.29 

29.62 

30.01 

30.36 

0.9277 

0.9294 

0.9021 

0.9409 

0.9412 

0.9252 

0.9113 

0.9274 

0.9215 

0.9329 

0.9260 

NCSR-GOC 

29.98 

31.03 

28.67 

34.31 

32.06 

30.04 

27.92 

30.18 

29.84 

30.29 

30.43 

0.9332 

0.9316 

0.9059 

0.9396 

0.9414 

0.9254 

0.9071 

0.9260 

0.9251 

0.9371 

0.9272 


Gaussian 

EISTA 

30.36 

29.36 

26.81 

31.50 

31.23 

29.47 

25.03 

29.65 

29.42 

29.36 

29.22 

0.9452 

0.9024 

0.8845 

0.8968 

0.9290 

0.9011 

0.8415 

0.9256 

0.9057 

0.9393 

0.9071 

ASDS (3 

29.83 

30.27 

27.29 

31.87 

32.93 

30.36 

27.05 

31.91 

28.95 

30.62 

30.11 

0.9126 

0.9064 

0.8637 

0.8978 

0.9576 

0.9058 

0.8881 

0.9491 

0.9039 

0.9304 

0.9115 

IDD-BM3D 

30.73 

31.68 

28.17 

34.08 

32.89 

31.45 

27.19 

31.66 

29.99 

31.40 

30.92 

0.9442 

0.9426 

0.9136 

0.9359 

0.9561 

0.9430 

0.8986 

0.9496 

0.9373 

0.9512 

0.9372 

NCSR [T] 

30.84 

31.37 

28.27 

33.69 

33.40 

31.17 

28.02 

32.23 

30.01 

31.62 

31.06 

0.9379 

0.9348 

0.9044 

0.9339 

0.9589 

0.9360 

0.9108 

0.9533 

0.9300 

0.9514 

0.9351 

NCSR-GOC 

31.32 

31.48 

28.44 

33.80 

33.45 

31.28 

27.45 

32.27 

30.27 

32.04 

31.18 

0.9486 

0.9413 

0.9153 

0.9375 

0.9594 

0.9429 

0.9014 

0.9554 

0.9389 

0.9587 

0.9399 


parameters for NCSR are selected as Dong et al. Q in order 
to maintain the consistency in the comparison. 


The PSNR measures of the reconstruction qualities are 
presented in Table |V] The results obtained with the image 
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Fig. 10: Test images for denoising: Lena, Monarch, Barbara, Boat, Cameraman (C. Man), 
Couple, Fingerprint (F Print), Hill, House, Man, Peppers, Straw. 

denoising algorithms SAPCA-BM3D LSSC ||44j; EPLL 
| [45| ; NCSR |[T| reported in Q for the same experiments 
are also given for the purpose of comparison. The overall 
performances of all algorithms are observed to be quite close, 
and the best average PSNR is given by SAPCA-BM3D at most 
noise levels. Nevertheless, the comparison between NCSR 
and NCSR-AGNN is more interesting, which shows that the 
proposed NCSR-AGNN algorithm yields a very similar per¬ 
formance to NCSR in denoising. A very slight improvement 
in average PSNR is obtained over NCSR at small noise levels, 
while this small advantage is lost at large noise levels. One can 
observe that the performance of NCSR-AGNN is better on the 
Monarch and Fingerprint images. This may be an indication 
that in such images with strong and oscillatory high-frequency 
textures, the patch manifold must have a particular geometry 
that is easier to identify under noise and the consideration of 
the geometry in assigning the similarities may help improve 
the denoising performance. 

VII. Conclusion 

In this paper, we have focused on the problem of selecting 
local subsets of training data samples that can be used for 
learning local models for image reconstruction. This study has 
been motivated by the observation that the Euclidean distance 
may not always be a good dissimilarity measure for comparing 
data samples lying on a manifold. We have proposed two 
methods for such data subset selection which take into account 
the geometry of the data assumed to lie on a manifold. 
Although the addressed problem has close links with manifold 
clustering, it differs by the fact that the goal here is not to 
obtain a partitioning of data, but instead select a local subset 
of training data that can be used for learning a good model 
for sparse reconstruction of a given input test sample. The 
performance of the methods has been demonstrated in a super¬ 
resolution application leading to a novel single-image super¬ 
resolution algorithm which outperforms reference methods, as 
well as in deblurring and denoising applications. 
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TABLE V: PSNR (in dB) results 
LSSC [4^; EPLL f^; NCSR fTl 


for the luminance components of denoised images for different denoising algorithms are reported in the following order: SAPCA-BM3D (43| ; 
; and NCSR with proposed AGNN. 



Methods 

Lena 

Monarch 

Barbara 

Boat 

C. Man 

Couple 

F. Print 

Hill 

House 

Man 

Peppers 

Straw 

Average 


SAPCA-BM3D 

38.86 

38.69 

38.38 

37.50 

38.54 

37.60 

36.67 

37.31 

40.13 

37.99 

38.30 

35.81 

37.98 


LSSC 

38.68 

38.53 

38.44 

37.34 

38.24 

37.41 

36.71 

37.16 

40.00 

37.84 

38.15 

35.92 

37.87 

a = 6 

EPLL 

38.52 

38.22 

37.56 

36.78 

38.04 

37.32 

36.41 

37.00 

39.04 

37.67 

37.93 

35.36 

37.49 


NCSR 

38.70 

38.49 

38.36 

37.35 

38.17 

37.44 

36.81 

37.17 

39.91 

37.78 

38.06 

35.87 

37.84 


NCSR-AGNN 

38.74 

38.62 

38.32 

37.34 

38.19 

37.40 

36.86 

37.15 

40.06 

37.78 

38.09 

35.82 

37.86 


SAPCA-BM3D 

36.07 

34.74 

35.07 

34.10 

34.52 

34.13 

32.65 

33.84 

37.06 

34.18 

34.94 

31.46 

34.40 


LSSC 

35.83 

34.48 

34.95 

33.99 

34.14 

33.96 

32.57 

33.68 

37.05 

34.03 

34.80 

31.39 

34.24 

a = 10 

EPLL 

35.56 

34.27 

33.59 

33.63 

33.94 

33.78 

32.13 

33.49 

35.81 

33.90 

34.51 

30.84 

33.79 


NCSR 

35.81 

34.57 

34.98 

33.90 

34.12 

33.94 

32.70 

33.69 

36.80 

33.96 

34.66 

31.50 

34.22 


NCSR-AGNN 

35.84 

34.66 

34.94 

33.87 

34.13 

33.90 

32.72 

33.66 

36.87 

33.95 

34.69 

31.46 

34.22 


SAPCA-BM3D 

34.43 

32.46 

33.27 

32.29 

32.31 

32.20 

30.46 

32.06 

35.31 

32.12 

33.01 

29.13 

32.42 


LSSC 

34.14 

32.15 

32.96 

32.17 

31.96 

32.06 

30.31 

31.89 

35.32 

31.98 

32.87 

28.95 

32.23 

a = lb 

EPLL 

33.85 

32.04 

31.33 

31.89 

31.73 

31.83 

29.83 

31.67 

34.21 

31.89 

32.56 

28.50 

31.78 


NCSR 

34.09 

32.34 

33.02 

32.03 

31.99 

31.95 

30.46 

31.86 

35.11 

31.89 

32.70 

29.13 

32.21 


NCSR-AGNN 

34.11 

32.37 

32.98 

32.01 

32.00 

31.94 

30.47 

31.84 

35.14 

31.88 

32.73 

29.14 

32.22 


SAPCA-BM3D 

33.20 

30.92 

31.97 

31.02 

30.86 

30.83 

28.97 

30.85 

34.03 

30.73 

31.61 

27.52 

31.04 


LSSC 

32.88 

30.58 

31.53 

30.87 

30.54 

30.70 

28.78 

30.71 

34.16 

30.61 

31.47 

27.36 

30.85 

o 

(M 

II 

b 

EPLL 

32.60 

30.48 

29.75 

30.63 

30.28 

30.47 

28.29 

30.47 

33.08 

30.53 

31.18 

26.93 

30.39 


NCSR 

32.92 

30.69 

31.72 

30.74 

30.48 

30.56 

28.99 

30.61 

33.97 

30.52 

31.26 

27.50 

30.83 


NCSR-AGNN 

32.89 

30.72 

31.70 

30.73 

30.50 

30.55 

29.01 

30.52 

33.98 

30.51 

31.28 

27.50 

30.82 


SAPCA-BM3D 

29.07 

26.28 
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o 

UT) 

II 

b 
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25.67 
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26.22 
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26.60 
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24.53 

26.86 

29.63 

26.60 

26.53 

22.48 

26.44 


NCSR-AGNN 

28.90 

25.69 

27.08 

26.57 

26.12 

26.19 

24.50 

26.80 

29.63 

26.59 

26.54 

22.46 

26.42 


SAPCA-BM3D 

25.37 

22.31 

23.05 

23.71 

22.91 

23.19 

21.07 

24.10 

25.20 

23.86 

23.05 

19.42 

23.10 


LSSC 

25.96 

21.82 

23.56 

23.94 

23.14 

23.34 

21.18 

24.30 

25.63 

24.00 

23.14 

19.50 

23.29 

fj = 100 

EPLL 

25.30 

22.04 

22.10 

23.78 

22.87 

23.34 

19.80 

24.37 

25.44 

23.96 

22.93 

18.95 

22.91 


NCSR 

25.66 

22.05 

23.30 

23.64 

22.89 

23.22 
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